Assignment 3: Methods of Moment Estimation

Assignment Objectives

• Master the fundamental concepts of point estimation and performance metrics

• Understand the theoretical foundation of the method of moments estimator (MME)

• Implement MME in R, incorporating numerical approximation methods

Log-logistic Distribution

The log-logistic distribution (also known as the Fisk distribution) is a continuous probability distribution that is particularly useful in contexts where data exhibit non-negative, skewed behavior and where the hazard rate is unimodal (increases to a peak and then decreases). It has been widely used in the areas such as survival analysis and reliability engineering, environmental science, economics, pharmacology, finance and risk management, etc.

For given shape parameter \(\beta\) and scale parameter \(\alpha\), the cumulative distribution function

\[ F(x) = \frac{1}{1+(x/\alpha)^{-\beta}} \]

As an exercise, you can derive the density in the following form

\[ f(x) = \frac{(\beta/\alpha)(x/\alpha)^{\beta-1}}{[1+(x/\alpha)^\beta]^2}, \ \ \text{ for } \ \ x > 0. \]

After some algebra, we can find the \(k\)th moment

\[ \mu_k = E[X^k] = \alpha^k B\left(1+\frac{k}{\beta}, 1 - \frac{k}{\beta} \right). \]

This assignment will focus on finding MME of parameters \(\alpha\) and \(\beta\) based on a real-world application data set.

Question 1: Derive the log-logistic density function

Given the CDF of the two-parameter log-logistic distribution

\[ F(x) = \frac{1}{1+(x/\alpha)^{-\beta}}. \] $$ \[\begin{aligned} F(x) &= \frac{1}{1 + (x/\alpha)^{-\beta}} \\[6pt] &= \left(1 + (x/\alpha)^{-\beta}\right)^{-1} \\[10pt] f(x) &= \frac{d}{dx}F(x) \\[6pt] &= -\left(1 + (x/\alpha)^{-\beta}\right)^{-2} \cdot \frac{d}{dx}\left(1 + (x/\alpha)^{-\beta}\right) \\[10pt] &= -\left(1 + (x/\alpha)^{-\beta}\right)^{-2} \cdot \left(0 + (-\beta)(x/\alpha)^{-\beta-1}\cdot \frac{1}{\alpha}\right) \\[10pt] &= \frac{\beta (x/\alpha)^{-(\beta+1)}} {\alpha \left(1 + (x/\alpha)^{-\beta}\right)^2} \\[10pt] &= \frac{\beta (x/\alpha)^{\beta-1}} {\alpha \left(1 + (x/\alpha)^{\beta}\right)^2} \\[12pt] &= \frac{(\beta/\alpha)(x/\alpha)^{\beta-1}}{[1+(x/\alpha)^\beta]^2}, \ \ \text{ for } \ \ x > 0. \end{aligned}\]

$$

Question 2: Distribution of Recovery Time from A Surgery

Time to recovery (in days) after a specific knee surgery procedure. This follows a typical log-logistic pattern in medical survival/recovery analysis:

8.23, 12.74, 14.83, 16.61, 18.16, 19.55, 20.80, 21.94, 23.00, 23.98, 24.89, 25.75, 26.56, 
27.34, 28.08, 28.79, 29.48, 30.15, 30.81, 31.45, 32.08, 32.70, 33.31, 33.92, 34.53, 35.13, 
35.73, 36.33, 36.93, 37.53, 38.14, 38.75, 39.37, 40.00, 40.64, 41.29, 41.95, 42.63, 43.33, 
44.05, 44.79, 45.56, 46.36, 47.20, 48.08, 49.02, 50.03, 51.12, 52.32, 53.65

Based on the above data to perform the following analysis.

1. Using method of moment estimation to estimate \(\alpha\) and \(\beta\), denoted by \(\hat{\alpha}\) and \(\hat{\beta}\), respectively.

Moment generating function for a log logistic distribution. \[ \mu_k = E[X^k] = \alpha^k B\left(1+\frac{k}{\beta},\, 1-\frac{k}{\beta}\right) \] The first two moments are \[ \mu_1 = E[X] = \alpha B\left(1+\frac{1}{\beta},\, 1-\frac{1}{\beta}\right) \] \[ \mu_2 = E[X^2] = \alpha^2 B\left(1+\frac{2}{\beta},\, 1-\frac{2}{\beta}\right) \] Beta Function Identity \[ B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} \] Simplifying both moments \[ \mu_1 = \alpha \Gamma\left(1+\frac{1}{\beta}\right) \Gamma\left(1-\frac{1}{\beta}\right) \] \[ \mu_2 = \alpha^2 \Gamma\left(1+\frac{2}{\beta}\right) \Gamma\left(1-\frac{2}{\beta}\right) \] The first two sample moments are \[ m_1 = \frac{1}{n}\sum_{i=1}^n X_i = 34.1922 \] \[ m_2 = \frac{1}{n}\sum_{i=1}^n X_i^2 = 1288.845 \] Substitute theoretical sample moments for moments since they’re equal \[ 34.1922 = \alpha \Gamma\left(1+\frac{1}{\beta}\right) \Gamma\left(1-\frac{1}{\beta}\right) \]

\[ 1288.845 = \alpha^2 \Gamma\left(1+\frac{2}{\beta}\right) \Gamma\left(1-\frac{2}{\beta}\right) \] Solve first equation for \(\alpha\) \[ \alpha = \frac{34.1922} {\Gamma\left(1+\frac{1}{\beta}\right) \Gamma\left(1-\frac{1}{\beta}\right)} \]

Plug \(\alpha\) into second equation \[ 1288.845 = \left( \frac{34.1922} {\Gamma\left(1+\frac{1}{\beta}\right) \Gamma\left(1-\frac{1}{\beta}\right)} \right)^2 \Gamma\left(1+\frac{2}{\beta}\right) \Gamma\left(1-\frac{2}{\beta}\right) \]

Divide both sides by \({m_1^2}\) (replacing theoretical sample moments with \({m_1}\) and \({m_2}\)) \[ \frac{m_2}{m_1^2} = \frac{ \Gamma\left(1+\frac{2}{\beta}\right) \Gamma\left(1-\frac{2}{\beta}\right) }{ \left[ \Gamma\left(1+\frac{1}{\beta}\right) \Gamma\left(1-\frac{1}{\beta}\right) \right]^2 } \]

Define \[ g(\beta) = \frac{ \Gamma\left(1+\frac{2}{\beta}\right) \Gamma\left(1-\frac{2}{\beta}\right) }{ \left[ \Gamma\left(1+\frac{1}{\beta}\right) \Gamma\left(1-\frac{1}{\beta}\right) \right]^2 } - \frac{m_2}{m_1^2} \] Solve \(g(\beta)=0\) to obtain \(\hat{\beta}\), and then substitute into \[ \hat{\alpha} = \frac{m_1} {\Gamma\left(1+\frac{1}{\hat{\beta}}\right) \Gamma\left(1-\frac{1}{\hat{\beta}}\right)} \]

Code to solve \(g(\beta)=0\) and obtain both \(\hat{\beta}\) and \(\hat{\alpha}\).

x <- c(8.23, 12.74, 14.83, 16.61, 18.16, 19.55, 20.80, 21.94, 23.00, 23.98, 24.89, 25.75, 26.56, 
27.34, 28.08, 28.79, 29.48, 30.15, 30.81, 31.45, 32.08, 32.70, 33.31, 33.92, 34.53, 35.13, 
35.73, 36.33, 36.93, 37.53, 38.14, 38.75, 39.37, 40.00, 40.64, 41.29, 41.95, 42.63, 43.33, 
44.05, 44.79, 45.56, 46.36, 47.20, 48.08, 49.02, 50.03, 51.12, 52.32, 53.65)

#1st and 2nd moments
m1 = mean(x)
m2 = mean(x^2)

#root finding equation
gb <- function(B){
  ( (gamma(1 + 2/B) * gamma(1 - 2/B)) / ((gamma(1 + 1/B) * gamma(1 - 1/B))^2) ) - (m2 / m1^2)
}

# since beta > 2
x_vals <- seq(2.01, 10, length.out = 200)
dframe <- data.frame(x = x_vals, y = gb(x_vals))

gk.plot <- ggplot(dframe, aes(x = x, y = y)) +
  geom_line(color = "steelblue", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +  # x-axis
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +  # y-axis
  labs(title = "The curve of function g(Beta)",
       x = "Beta", 
       y = "g(Beta)") +
    theme(plot.title = element_text(hjust = 0.5),
        plot.margin = margin(t = 35, r = 20, b = 30, l = 30, unit = "pt"))
ggplotly(gk.plot)

#find the root (beta value when equation = 0)
B <- uniroot(gb, interval = c(5, 7))$root

#use the beta to find the alpha
alpha <- m1/ (gamma(1+1/B) * gamma(1-1/B))
pander(cbind(beta = B, alpha = alpha))

beta	alpha
6.006	32.65

2. Since the moment estimates \(\hat{\alpha}\) and \(\hat{\beta}\) are random, construct bootstrap sampling distributions for each. To visualize these distributions, plot separate bootstrap histograms for \(\hat{\alpha}\) and \(\hat{\beta}\). hen, overlay a smooth density curve on each histogram using Gaussian kernel density estimation. Finally, describe the patterns of these density curves.

bootstrap <- function (data, statistics, B){
  n <- length(data)
  stats <- matrix(nrow = B, ncol = 2)
  for (i in 1:B){
    new_data <- sample(data, n, replace = TRUE)
    stats[i,] <- statistics(new_data)
  }
  return(stats)
}
  
statistics <- function (data){
  m1 = mean(data)
  m2 = mean(data^2)
  
  gb <- function(B){
  ( (gamma(1 + 2/B) * gamma(1 - 2/B)) / ((gamma(1 + 1/B) * gamma(1 - 1/B))^2) ) - (m2 / m1^2)}
  
  x_vals <- seq(2.01, 200, length.out = 100)
  dframe <- data.frame(x = x_vals, y = gb(x_vals))
  
  B <- uniroot(gb, interval = c(3, 15))$root
  alpha <- m1/ (gamma(1+1/B) * gamma(1-1/B))
  
  c(B, alpha)
}

stats <- bootstrap(x, statistics, 5000)

beta_boot  <- stats[, 1]
alpha_boot <- stats[, 2]

hist(beta_boot,
     breaks = 20,
     probability = TRUE,
     xlab = expression(hat(beta)),
     main = "Bootstrap Sampling Distribution\n of Beta Hat",
     cex.main = 0.9,
     col = "lightgray",
     border = "white")

lines(density(beta_boot), col = "navy", lwd = 2)

The bootstrap sampling distribution of \(\hat{\beta}\) is roughly bell-shaped with a little bit of skewness to the right. It’s centered close to the original MoM estimate, which suggests the estimator is not heavily biased. The spread of the histogram shows how much \(\hat{\beta}\) would vary if we collected new samples from the population. Since the distribution isn’t extremely wide or irregular, it suggests that the MoM estimator for \({\beta}\) is fairly stable for this sample.

hist(alpha_boot,
     breaks = 20,
     probability = TRUE,
     xlab = expression(hat(alpha)),
     main = "Bootstrap Sampling Distribution\n of Alpha Hat",
     cex.main = 0.9,
     col = "lightgray",
     border = "white")

lines(density(alpha_boot), col = "darkred", lwd = 2)

The bootstrap sampling distribution of \(\hat{\alpha}\) is approximately symmetric and bell-shaped. It’s centered pretty close to the original MoM estimate, which suggests the estimator is not heavily biased. The spread of the histogram shows how much \(\hat{\alpha}\) would vary if we collected new samples from the population. Since the distribution isn’t extremely wide or irregular, it suggests that the MoM estimator for \({\alpha}\) is fairly stable for this sample.